/* Copyright (c) CERN
 * * Artistic License 2.0
 * */

import java.util.Random;
import java.lang.Math;

/**
 * Gamma distribution.
 * This is the only one that does not implement the MoM estimation.
 * It might be the reason why refinement does not seem to work.
 *
 * @author Xavier Grehant
 */
public class Gamma extends DistroUtils implements Distro {

/**
 * Shape
 */
   double k;

/**
 * Scale
 */
   double theta;

/**
 * A uniform random number generator, initialized at contruction.
 */
   Random rand;

/**
 * Used in the random generation with the acceptance-rejection method.
 */
   double delta;

/**
 * Used in the random generation with the acceptance-rejection method.
 */
   double v0;

/**
 * Constructor taking distribution parameters.
 * @param k shape.
 * @param theta scale.
 */
   public Gamma(double k, double theta) {
      this.k = k;
      this.theta = theta;
      setCrap(k, theta);
      rand = new Random();
   }

/**
 * This method computes the variables needed for the random number generation.
 * @return
 * @param k shape.
 * @param theta scale.
 */
   void setCrap(double k, double theta) {
      delta = k - (double) (int) k;
      v0 = Math.E / (Math.E + delta);
   }

/**
 * Constructor taking samples.
 * This constructor fits the distribution.
 * See http://en.wikipedia.org/wiki/Gamma_distribution
 * @param array an array of integer samples.
 */
   public Gamma(int[] array) {
      rand = new Random();
      array = discard(array, .01, 10);
      long cumul = 0, logcumul = 0;
      for (int i = 0; i < array.length; i++) {
          cumul += array[i];
          logcumul += Math.log(array[i]);
      }

      double s = Math.log((double) cumul/(double) array.length) - logcumul / (double) array.length;
  
      k = (3d - s + Math.sqrt(24d * s + Math.pow(s - 3d, 2d))) / (12d * s);
      for (int j = 0; j < 11; j++) {
         System.out.println(" k = " + k);
         if (Double.isNaN(k)) break;
         k = k - ((Math.log(k) - digamma(k) - s) / ((1d / k) - trigamma(k)));
      }
  
      theta = (double) cumul / (k * (double) array.length);
      System.out.println("k : " + k);
      System.out.println("s : " + s);
      System.out.println("theta : " + theta);
      setCrap(k, theta);
   }

/**
 * This factory method creates a Gamma distribution instance from a sample set.
 * It estimates the biased fit (see PDPTA'08 paper from Grehant and Jarp).
 * @return a Gamma distribution
 * @param values a sample set.
 * @param minvalue the minimum value of interest.
 * @param refinements the number of times the fit is refined.
 */
   static Gamma biased(int[] values, int minvalue, int refinements) {
      Gamma gamma = new Gamma(values);
      for (int refined = 0; refined < refinements; refined++) {
         values = gamma.bias(values, minvalue);
         gamma = new Gamma(values);
      }
      return gamma;
   }

/**
 * Random number generator according to the distribution.
 * @return a new random double according to the distribution.
 */
   public double nextDouble() {
      return theta * (xi() - sumlog());
   }

/**
 * Sum of sample logarithms
 * @return the sum of 'k' logarithms of random samples.
 */
   double sumlog() {
      double sum = 0;
      for (int i = 0; i < (int) k; i++)
         sum += Math.log(rand.nextDouble());
      return sum;
   }

/**
 * Variable used in the random number generation.
 * @return xi from uniform numbers.
 */
   double xi() {
      return xi(rand.nextDouble(), rand.nextDouble(), rand.nextDouble());
   }

/**
 * Acception-rejection method.
 * @return xi from given seeds.
 * @param v1 value at step n in the algorithm
 * @param v2 value at step n+1 in the algorithm
 * @param v3 value at step n+3 in the algorithm
 */
   double xi(double v1, double v2, double v3) {
      double xi, eta;
      if (v1 <= v0) {
         xi = Math.pow(v2, 1d / delta);
         eta = v3 * Math.pow(xi, delta - 1d);
      }
      else {
         xi = 1d - Math.log(v2);
         eta = v3 * Math.exp(-xi);
      }
      if (eta > Math.pow(xi, delta - 1d) * Math.exp(-xi))
         return xi(v2, v3, rand.nextDouble());
      return xi;
   }

/**
 * DiGamma function.
 * @return digamma function.
 * @param k shape.
 */
   static double digamma(double k) {
      if (k >= 8)
         return Math.log(k)
            - (1d
                  + ((1d
                        - (((1d / 10d)
                           - (1d
                              / (21d
                                 * Math.pow(k, 2d))))
                     / Math.pow(k, 2)))
               / (6 * k)))
         / (2d * k);
      else {
         return digamma(k + 1d) - (1d / k);
      }
   }

/**
 * TriGamma function.
 * @return trigamma function.
 * @param k shape.
 */
   static double trigamma(double k) {
      if (k >= 8)
         return
      (1d
            + (1d
                  + ((1d
                        - (((1d / 5d)
                           - (1d
                              / (7d
                                 * Math.pow(k, 2d))))
                     / Math.pow(k, 2)))
               / (3 * k)))
         / (2d * k))
   / k;
      else
         return trigamma(k + 1d) + (1d / Math.pow(k, 2));
   }
}
